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Abstract.  The  Taylor-Couette  problem  for  a  rarefied  gas  is  studied  numerically  by  the  direct  simulation  Monte  Carlo  method. 
The  gas  is  supposed  to  be  contained  in  an  annular  domain,  bounded  by  two  coaxial  rotating  cylinders  and  top  and  bottom 
(specularly  reflecting)  boundaries,  and  the  flow  is  assumed  to  be  steady  and  axisymmetric.  Special  attention  is  focused  on  the 
effect  of  rotation  of  the  outer  cylinder  on  the  type  of  the  induced  Taylor- vortex  flow.  It  is  shown  that  different  types  of  flow 
can  coexist  stably  in  a  wide  range  of  speeds  of  rotation  of  the  inner  and  outer  cylinders  unless  the  outer  cylinder  is  rotating 
fast  in  the  opposite  direction  to  the  inner. 


INTRODUCTION 

The  Taylor-Couette  problem,  the  instability  of  the  cylindrical  Couette  flow  and  the  occurrence  of  the  secondary 
flow  with  vortices  (Taylor  vortex  flow)  in  a  fluid  between  two  coaxial  rotating  cylinders,  is  a  classical  problem  in 
fluid  dynamics  (e.g..  Ref.  [1,  2,  3]).  However,  the  study  of  this  problem  for  a  rarefied  gas  started  much  later,  and 
several  numerical  studies  with  different  emphases  [4,  5,  6,  7,  8]  have  been  carried  out  by  using  the  direct  simulation 
Monte  Carlo  (DSMC)  method  [9].  For  instance,  the  effect  of  the  temperature  difference  between  the  two  cylinders  is 
investigated  in  detail  in  Ref.  [8].  Lately,  interest  has  been  expanded  to  more  complicated  topics,  such  as  the  bifurcation 
of  flows  when  evaporation  and  condensation  take  place  on  the  cylinders  [10]  and  the  ghost  effect  when  the  rotation 
speeds  of  the  cylinders  as  well  as  the  Knudsen  number  are  infinitesimal  [11]. 

In  the  previous  DSMC  studies,  the  case  where  only  the  inner  cylinder  is  rotating  is  mainly  discussed.  In  the 
present  paper,  we  investigate  the  same  problem  again  by  the  DSMC  method,  restricting  ourselves  to  the  steady  and 
axisymmetric  flows.  But  our  attention  is  focused  on  systematic  understanding  of  the  effect  of  rotation  of  the  outer 
cylinder  on  the  type  of  induced  Taylor  vortex  flow  and  the  parameter  range  in  which  different  types  of  flow  may  exist. 

PROBLEM  AND  ASSUMPTIONS 

Let  us  consider  a  rarefied  gas  between  two  rotating  coaxial  circular  cylinders  with  a  finite  length.  To  describe  the 
problem,  we  use  the  cylindrical  coordinate  system  (r,  6,  z)  with  the  z  axis  taken  along  the  common  axis  of  the  cylinders. 
The  radius  of  the  inner  cylinder  is  L/,  that  of  the  outer  cylinder  is  Lu,  and  the  bottom  and  top  ends  of  the  cylinders 
are  covered  with  plates  located  at  z  =  0  and  Lz,  respectively.  Thus,  we  consider  an  annular  domain  Lj  <  r  <  Lu, 
0  <  6  <  2 n,  and  0  <  z  <  Lz.  The  inner  and  outer  cylinders  are  rotating  around  the  z-axis  at  surface  velocities  V;  and 
Vu  in  the  6  direction,  respectively,  and  both  cylinders  are  kept  at  a  common  temperature  Tq.  We  investigate  the  steady 
behavior  of  the  gas  numerically  on  the  basis  of  kinetic  theory  under  the  following  assumptions:  (i)  the  flow  field  is 
axisymmetric;  (ii)  the  gas  molecules  are  hard  spheres;  (iff)  the  gas  molecules  undergo  diffuse  reflection  on  the  surface 
of  the  cylinders  (r  =  Lj  and  Lu)  and  specular  reflection  on  the  bottom  and  top  boundaries  (z  =  0  and  L-).  In  the  present 
study,  special  attention  is  focused  on  the  effect  of  the  rotation  of  the  outer  cylinder  on  the  flow  patterns. 

Since  the  basic  equations  are  essentially  the  same  as  those  of  Ref.  [8]  and  the  space  is  limited  in  the  present 
proceedings,  we  omit  them  for  brevity.  We  only  note  that,  by  introducing  appropriate  dimensionless  variables,  one  finds 
that  the  present  problem  is  characterized  by  the  following  five  dimensionless  parameters:  L///L/,  L-/L/,  V// {2RTf) 1  1 , 
Vu/ (2RTo) 1 Y  and  Kn.  Here,  Kn=  £q/Lj  is  the  Knudsen  number,  where  £q  is  the  mean  free  path  of  the  gas  molecules 
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FIGURE  1.  Map  of  the  flow  patterns  in  the  (Vji / (IRTq)1  /2 ,  Vj  / ( 2RTq )1/2)  plane  for  L///L/  =  2,  Lz/Li  =  1,  and  Kn  =  0.01.  The 
Maxwellian  distribution  corresponding  to  the  equilibrium  state  at  rest  with  temperature  Tq  and  density  pav  is  used  as  the  initial 
condition.  The  symbols  o,  •,  and  T  indicate  that  a  Couette  flow,  a  single-vortex  flow,  and  a  double-vortex  flow,  respectively,  are 
obtaind  at  the  corresponding  points. 


in  the  equilibrium  state  at  rest  with  temperature  Tq  and  density  pm,  (pav:  average  density  of  the  gas  over  the  annular 
domain),  and  R  is  the  gas  constant  per  unit  mass. 

As  the  solution  method,  we  use  the  DSMC  method  by  Bird  [9].  As  in  the  usual  DSMC  computation,  we  obtain 
the  steady  flow  field  as  the  long-time  limit,  pursuing  the  long-time  behavior  of  the  solution  of  the  time-dependent 
boundary-value  problem  with  an  appropriately  chosen  initial  condition. 

RESULTS  AND  DISCUSSIONS 

Since  the  DSMC  method  is  explained  in  many  places,  we  omit  its  description  for  short  and  show  some  of  the  results 
of  computation.  More  specifically,  we  restrict  ourselves  to  the  case  where  L///L/  =  2,  L-/ Lj  =  1,  and  Kn  =  0.01  and 
investigate  the  effect  of  rotation  speeds  V/ /  (2RTq)  1  2  and  Vn/ ( IRTq )  of  the  cylinders  on  the  flow  patterns. 

In  the  actual  computation,  the  square  cross  section  (L/  <  r  <  Ljj,  0  <z<L-  with L// / Lj  =  2  and Lz/Lj  =  1)  is  divided 
into  120  x  120  (System  I)  or  240  x  240  (System  II)  square  cells  of  a  uniform  size;  the  total  number  of  simulation 
particles  is  144  x  104  (System  I)  or  576  x  104  (System  II);  and  the  time  step  At  is  At  —  5  x  10  4L/(2RTq)  1/2  = 
5  x  10_4(2/y/7t)toKn_1 ,  where  fo  is  the  mean  free  time  corresponding  to  £q .  Except  Fig.  7,  the  results  presented  in 
this  paper  are  based  on  System  I. 

First  we  perform  computation,  taking  the  Maxwellian  distribution  corresponding  to  the  equilibrium  state  at  rest 
with  temperature  Tq  and  density  pav  as  the  initial  condition  at  time  t  =  0.  This  means  that  the  cylinders  start  to  rotate 
impulsively  at  t  =  0.  The  computation  for  various  values  of  V/ / (2RTq)  '/2  and  Vn / (2RTq)]  -  gives  three  different  types 
of  steady  flow  pattern:  an  axially  uniform  flow  (Couette  flow);  a  Taylor  vortex  flow  with  a  single  vortex  (single- vortex 
flow);  and  a  Taylor  vortex  flow  with  double  vortices  arranged  in  the  axial  direction  (double-vortex  flow).  A  map  of  the 
flow  patterns  in  the  (V u  j  (2RTq)  1  /2.  V/  /  (2R'Tq)  l/2)  plane  is  shown  in  Fig.  1.  A  Couette  flow  is  obtained  at  the  points 
with  o,  a  single-vortex  flow  at  the  points  with  •,  and  a  double-vortex  flow  at  the  points  with  ▼.  The  negative  Vn  means 
that  the  outer  cylinder  is  rotating  in  the  direction  opposite  to  the  inner. 

For  a  fixed  Vn / (2RTq)  '/2,  a  Couette  flow  is  obtained  for  Vi  /  (2RTq)  1  2  smaller  than  a  critical  value.  The  critical  value 
is  smallest  at  around  Vn / (2RTq)  12  =  —0.1  and  0  and  increases  rapidly  with  the  increase  of  V// / (2RTq ) 1  2.  When  the 
rotation  speed  of  the  outer  cylinder  is  slow,  a  single-vortex  flow  is  obtained  in  a  rather  wide  range  of  Vj/ (2RTq) 
However,  this  range  shrinks  and  a  double-vortex  solution  becomes  dominant  as  the  rotation  speed  of  the  outer  cylinder 
increases,  in  particular,  in  the  opposite  direction. 

Figure  1  also  shows  an  isolated  double-vortex  flow  at  V/ / (2RTq)  1  =  0.9  in  the  middle  of  the  range  of  a  single¬ 
vortex  flow  for  Vii/(2RTq)1'2  =  —0.1.  This  fact  suggests  the  possibility  of  the  coexistence  of  a  single-  and  a  double¬ 
vortex  flow  around  this  parameter  range.  Our  next  interest  is  to  investigate  whether  there  is  a  parameter  range  in  which 
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FIGURE  2.  Existence  range  of  the  steady  solutions  in  the  (Vji  / {2RTg)^/2 ,  Vj  /  (2RTg)1/2)  plane  for  Ljj  /Lj  =  2,  Lz/Lj  =  1,  and 
Kn  =  0.01.  The  symbols  o,  •,  and  Y  indicate  that  a  Couette  flow,  a  single-vortex  flow,  and  a  double-vortex  flow,  respectively,  are 
obtained  at  the  corresponding  points.  Here,  the  symbols  stand  for  the  data  precisely  at  Vjj / (2RTg)1/2  =  —0.6,  —0.3,  —0.1,  0,  0.1, 
0.3,  and  0.6,  but  »’s  (v’s)  are  shifted  leftward  (rightward)  slightly  for  legibility. 


different  types  of  solution  coexist  stably. 

For  this  purpose,  we  carry  out  computation  with  initial  conditions  different  from  the  Maxwellian  used  to  obtain  the 
data  in  Fig.  1.  For  example,  for  Vji  /  (2RTg)1^2  =  —0.1,  we  take  the  isolated  double-vortex  solution  at  Vi  /  (2RTg) l/'2  = 
0.9  as  the  initial  condition  for  a  (slightly)  smaller  value  of  V//(2R'i}))i/2,  say  V/  /  (2/(71)) 1/2  =  0.8.  If  this  computation 
leads  to  a  double-vortex  flow,  then  we  use  it  as  the  initial  condition  for  a  smaller  value  of  Vj/(2RTg)1/2,  say 
VI/(2RT0)1'2  =  0.7.  We  continue  this  process  until  the  computation  leads  to  a  solution  of  a  different  type  (i.e.,  a 
single-vortex  flow  or  a  Couette  flow)  or  a  fluctuating  state  that  cannot  be  regarded  as  a  steady  solution.  We  take  a 
similar  procedure  for  V/ / (2RTq)  1  2  larger  than  0.9.  In  this  way  one  can  try  to  extend  the  range  of  a  stable  double¬ 
vortex  flow  for  V/i  j (2RTq)  1  2  =  —0.1  as  wide  as  possible.  We  carry  out  the  same  procedure  to  extend  the  range  of  a 
stable  single-vortex  solution,  starting  from  that  at  the  edge  of  the  range  in  Fig.  1  [i.e.,  Vj/ (2RTg) 1  2  =  0.45  and  1.1],  as 
well  as  the  range  of  a  stable  Couette  flow  to  larger  values  of  Vj /  (2 RTq)1/2.  Essentially  the  same  procedure  is  carried  out 
for  other  values  of  Vn/(2RTg)i  2.  The  result  obtained  by  the  above  procedure  is  shown  in  Fig.  2.  The  symbols,  which 
have  the  same  meaning  as  in  Fig.  1,  indicate  the  results  at  Vji/ (2RTo) 1  </2  =  —0.6,  —0.3,  —0.1,  0,  0.1,  0.3,  and  0.6  as 
in  Fig.  1,  but  »'s  (y’s)  are  shifted  leftward  (rightward)  slightly  for  legibility.  In  the  process  of  extending  the  Couette  - 
flow  region,  we  obtained  a  fluctuating  state  that  cannot  be  regarded  as  a  steady  solution  at  some  points  at  the  edge  of 
the  region,  but  they  are  not  shown  in  the  figure.  There  is  a  significant  overlap  between  the  region  of  a  single-vortex 
flow  and  that  of  a  double-vortex  flow  when  V7// / (2/(71)) 1  /2  >  —0.1.  For  V// / (2RTg ) 1  /2  =  —0.3,  a  single-vortex  flow  is 
limited  to  a  rather  narrow  range  of  intermediate  Vi /  (2RTq)1'2 ,  and  its  overlap  with  the  range  of  a  double- vortex  flow  is 
small.  There  is  no  region  of  a  single -vortex  flow  for  Vii/(2RTq)1/2  =  —0.6  [in  this  case,  we  carried  out  computations 
for  Vj / (2RTq)1/2  =  0.75,  8,  9,  and  1,  using  the  single-vortex  flow  at  Vj  /  {2RTg)ll2  =  0.75  and  Vjj / (2RTq)1/2  =  —0.3 
as  the  initial  condition  and  obtained  only  steady  double-vortex  flows]. 

In  Figs.  3-5,  typical  flow  patterns  of  a  single-  and  a  double-vortex  flow  are  shown  for  three  points  in  Fig.  3  where  the 
two  types  of  flow  coexist,  i.e.,  (V)//(2/?7b)1/2,  V//(2/?7’o)1/2)  =  (0,  1.5)  in  Fig.  3,  (0.6,  3)  in  Fig.  4,  and  (—0.3,  0.75) 
in  Fig.  5.  The  panels  (a)  correspond  to  the  single-vortex  flow,  and  the  panels  (b)  to  the  double -vortex  flow.  The  arrow 
indicates  the  flow  velocity  vector  (vr,  v-)  in  the  (r,  z)  plane  at  its  starting  point,  where  vr  and  v-  are,  respectively,  the 
r  and  z  components  of  the  flow  velocity,  and  the  scale  of  (v2  +  v2)1/2/(2/?7’o)1/2  =  0.2  or  0.1  is  shown  in  each  figure. 
The  contour  line  vg  /  (2RTo) 1  2  =  const,  where  vg  is  the  6  component  of  the  flow  velocity,  is  also  shown  in  the  figures. 
The  obtained  double -vortex  flows  are  all  symmetric  with  respect  to  the  plane  z  =  L-/2  and  have  a  clockwise  vortex  in 
the  upper  half  (Lz/2  <  z<  Lz)  and  a  counterclockwise  vortex  in  the  lower  half  (0  <  z  <  Lz/2).  In  Fig.  5(a),  the  center 
of  the  vortex  is  as  low  as  that  of  the  lower  vortex  in  Fig.  5(b),  and  there  is  almost  no  flow  in  the  upper  part. 

A  convenient  physical  quantity  that  describes  the  difference  of  the  three  types  of  solution  quantitatively  would 
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FIGURE  3.  Flow  field  of  single-  and  double-vortex  flows  for  Vj / ( 2RTq )'/2  =  1.5  and  Vn/ (2/?7(j) 1  “  =  0  (L///L/  =  2,  Lz/Li  =  1, 
Kn  =  0.01 ).  (a)  single-vortex  flow,  (b)  double- vortex  flow.  The  arrow  indicates  the  flow  velocity  (vr,  vz )  in  the  (r,  z)  plane  at  its 
starting  point,  and  the  scale  of  (v2  +vz)1^2/(2RTq)1^2  =  0.2  is  shown  in  each  figure.  The  contour  lines  vg / (2RTg)1/2  =  const  are 
also  shown  in  the  figure. 


FIGURE  4.  Flow  field  of  single-  and  double-vortex  flows  for  Vj / (2 KTy)1/2  =  3  and  Vn  j ( 2RT§  j1/2  =  0.6  (L///L/  =  2,  Lz/Li  =  1, 
Kn  =  0.01 ).  (a)  single-vortex  flow,  (b)  double-vortex  flow.  See  the  caption  of  Fig.  3. 


FIGURE  5.  Flow  field  of  single-  and  double-vortex  flows  for  Vj/(2RTq )1/2  =  0.75  and  V///(27?7())1/2  =  —0.3  (L///L/  =  2, 
Lz/Lj  =  1,  Kn  =  0.01).  (a)  single-vortex  flow,  (b)  double-vortex  flow.  See  the  caption  of  Fig.  3.  Here,  the  scale  of  (v2  + 
V2)1/2/(2Rr0)1/2 

=  0.1  is  shown  in  each  figure. 
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FIGURE  6.  F2  versus  Vj/(2RTo )1/2  for  typical  values  of  V/j / (2RTy)1/2  for  L///Z./  =  2,  Lz/Lj  =  1,  and  Kn  =  0.01.  (a) 
Vji / (2RTq)1/2  =  0,  (b)  —0.1,  (c)  —0.3,  and  (d)  —0.6.  The  symbols  o,  •  ,  and  Y  indicate  a  Couette  flow,  a  single-vortex  flow, 
and  a  double-vortex  flow,  respectively. 


be  the  circulation  of  the  flow  velocity  vector  (tv,  v:)  in  the  (r.  z)  plane  [8,  10].  We  introduce  the  dimensionless 
circulation  Ti  and  r2  defined  by  F„,  =  am(2RTo)~  1//2/o27r[vc(fm,  zm)cos</>  —  vr(rm,  zm)  sin0]d0  (m  =  1,2),  where 
rm  =  L/(1.5  +  am cos <j> )  and  zm  =  Li(am  +  amsin(f))  with  a\  =  0.3,  a.\  =  0.5  and  a2  =  0.23,  a2  =  0.25;  that  is,  Ti 
(r2)  is  the  circulation  along  the  circle  of  radius  0.3 L/  (0.23 L/)  centered  at  (r,  z)  =  (1.5,  0.5 )Lj  [(r,  z)  =  (1.5,  0.25)L/]. 
Note  that  Tj  =  T2  =  0  for  a  Couette  flow,  and  F|  —  0  for  a  double-vortex  flow  that  is  symmetric  with  respect  to  the 
middle  plane  z  =  Lz/ 2.  Therefore,  to  distinguish  the  three  types  of  flow  at  the  same  time,  it  is  more  convenient  to  use 
r2.  Figure  6  shows  F 2  versus  V/ / ( 2/f7o ) 1  /  2  for  Vu / (2R1})) 1  /2  =  0,  —0.1,  —0.3,  and  —0.6.  Here,  the  values  of  F2  are 
calculated  using  the  averages  of  the  solution  over  5  x  105  time  steps  after  the  steady  state  is  established.  From  Fig.  6, 
it  is  likely  that  both  branches,  the  branch  of  a  single-vortex  flow  and  that  of  a  double-vortex  flow,  bifurcate  from  the 
Couette  flow.  However,  the  double-vortex  flow  may  be  unstable  for  the  values  of  F//(2/F/ii ) 1  2  near  the  bifurcation 
point  except  the  case  of  high-speed  opposite  rotation,  F///(2/F/o) 1  2  =  —0.6.  In  contrast,  the  single- vortex  flow  may 
lose  stability  rapidly  from  the  large  values  of  V//(2i?Fo)1/2  when  the  rotation  speed  of  the  outer  cylinder  increases  in 
the  opposite  direction,  and  it  may  become  unstable  for  all  Vj/ (2RTo) 1  /2  at  V//  j (2RTo) 1  /2  =  —0.6. 

Figures  2  and  6(b)  show  that  both  of  a  Couette  flow  and  a  single-vortex  flow  exist  stably  at  V//(27?Fo)1/2  =  0.4  in 
the  case  of  Vu/{2RTq)1/2  =  —0.1.  Such  an  overlap  between  the  range  of  a  Couette  flow  and  a  single-vortex  flow  is 
not  observed  for  other  V// / [2RT{) ) 1  2.  Therefore,  we  have  to  verify  whether  there  is  such  a  small  overlap  or  not  more 
carefully.  As  pointed  out  with  a  concrete  example  in  [10],  it  is  generally  very  difficult  to  obtain  reliable  results  for 
sensitive  problems,  such  as  the  bifurcation  of  flows  and  the  onset  of  instability,  by  the  DSMC  method  unless  one  uses 
very  fine  cells  and  an  extremely  large  number  of  simulation  particles.  Thus,  we  investigate  more  closely  the  behavior  of 
the  solution  in  the  vicinity  of  the  point  where  a  Couette  flow  and  a  single-vortex  flow  coexist,  i.e.,  V/ /(IRTo ) 1  2  =  0.4 
and  V/i / (2RTq)  ]  2  =  —0.1,  using  a  larger  computational  system  (System  II).  The  result  of  the  computation  is  shown  in 
Fig.  7(a),  where  T]  rather  than  F2  is  shown  because  the  difference  between  a  Couette  flow  and  a  single -vortex  flow  is 
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FIGURE  7.  Results  obtained  with  System  II  for  Vii/(2RTq j1/2  —  —0.1  (L///L/  =  2,  Lz/Lj  =  1,  Kn  =  0.01).  (a)  F)  vs 
Vi/(2RTqY/2.  The  symbols  o  and  •  indicate  a  Couette  flow  and  a  single-vortex  flow,  respectively,  (b)  Time  evolution  of  Fj 
for  F//(2Rr0)1/2  =  0.39.  The  symbol  ■  indicates  the  data  with  the  initial  condition  of  the  single-vortex  solution  obtained  for 
Vj / {2RTqY  '  -  =  0.40,  and  □  those  with  the  initial  condition  of  the  artificial  Couette  flow. 


clearer  for  T i .  The  single-vortex  solution  at  V/  /  (2RTq)  1  /  2  =  0.43  was  obtained  by  the  impulsive  rotation  as  in  Fig.  1 . 
Then  we  extended  the  branch  of  a  single -vortex  flow  down  to  Vj / (2RTo)  P2  =  0.40  by  the  same  way  as  the  construction 
of  Fig.  2.  However,  the  computation  for  Vi  /  (2RTq) 1  /2  =  0.39  using  the  single-vortex  solution  at  Vj  /  (2RTo) 1  '2  =  0.40  as 
the  initial  condition  did  not  converge  to  any  clear  steady  solution.  The  time  evolution  of  the  circulation  T i  in  this  case 
is  shown  by  the  symbol  ■  in  Fig.  7(b).  On  the  other  hand,  the  Couette  flows  at  Vj/ (2RTq) 1/  2  =  0.3,  0.35,  and  0.38  were 
obtained  from  special  initial  conditions  explained  in  the  following.  One  can  always  obtain  the  axially  uniform  Couette 
flow  by  imposing  the  uniformity  in  the  z  direction.  By  using  such  a  solution  in  each  layer  of  the  cells  piling  up  in 
the  z  direction,  we  can  generate  an  artificial  Couette  flow  without  constraint  of  axial  uniformity.  Such  flows  generated 
for  Vi/(2RTo)1'2  =  0.3,  0.35,  and  0.38  were  used  as  the  initial  conditions  for  corresponding  cases.  It  is  likely  that 
one  can  obtain  a  stable  Couette  flow  more  easily  from  this  type  of  initial  condition.  However,  such  a  computation  at 
Vi l (2RT()) 1  2  =  0.39  did  not  lead  to  any  clear  steady  flow.  The  time  evolution  of  F|  in  this  case  is  shown  by  the  symbol 
□  in  Fig.  7(b).  In  Fig.  7(b),  the  Fi  at  t/At  =  n  is  calculated  with  vr  and  v-  averaged  over  t/At  =  (n  —  5  x  104  +  1)  ~  n. 
By  the  way,  in  all  the  flows  that  we  judged  as  Couette  flows,  T]  and  F 2  at  large  time  exhibit  rather  high-frequency 
oscillations  around  zero  as  time  goes  on  .  From  Fig.  7,  we  may  conclude  that  there  is  no  overlap  between  the  range  of 
a  Couette  flow  and  that  of  a  single-vortex  flow  and  that  the  overlap  observed  in  Fig.  6(b)  is  a  false  behavior  caused  by 
the  insufficiency  of  the  computational  system  (System  I). 
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